subroutine build_exchange_lists(ispecies)
  use logs
  use my_mpi
  use parameter_module
  use camera_module, only : x0,y0,z0,L
  use particles    , only : NumPart_ByPE, readCoordinates, readSmoothingLength
  use grid
  use exchange
  implicit none
  integer*4 :: ispecies
  integer*4 :: i, iplane, ierr
  real*4, dimension(3) :: delta_coord,coord_min

  ! Allocate and initialise arrays for particle import and export
  call allocate_exchange_counters()
  NumRetain         = 0
  NumImportPE(:)    = 0
  NumExportPE(:)    = 0
  ExchangeMatrix(:) = 0
  ImportOffset(1)   = 0
  ExportOffset(1)   = 0

  ! Also allocate array for particle destination
  call allocate_destination_array()
  DestinationPE(:) = -1

  ! Construct export list 
  do i=1,NumPart_ByPE(MyPE)
     
     if(readCoordinates(1,i)+readSmoothingLength(i) < GridExtent(1,1)) cycle
     if(readCoordinates(1,i)-readSmoothingLength(i) > GridExtent(2,1)) cycle
     if(readCoordinates(2,i)+readSmoothingLength(i) < GridExtent(1,2)) cycle
     if(readCoordinates(2,i)-readSmoothingLength(i) > GridExtent(2,2)) cycle
     if(readCoordinates(3,i)+readSmoothingLength(i) < GridExtent(1,3)) cycle
     if(readCoordinates(3,i)-readSmoothingLength(i) > GridExtent(2,3)) cycle

#ifdef CUBE
     iplane = int((readCoordinates(3,i) - GridExtent(1,3)) / CellSize(3)) + 1
     if(iplane < 1 .or. iplane > params%nz) cycle
#else
     iplane = int((readCoordinates(3,i) - GridExtent(1,3)) / CellSize(2)) + 1
     if(iplane < 1 .or. iplane > params%nz) cycle
#endif

     if(iplane <= NZMaster)then
        ! Master process
        DestinationPE(i) = 1 
     else
        ! Slave processes
        iplane  = iplane - NZMaster - 1
        DestinationPE(i) = int(iplane/NZLocal) + 2
     endif
     
     if(DestinationPE(i) < 1 .or. DestinationPE(i) > NumPEs) cycle

     if(DestinationPE(i) == (MyPE+1))then
        NumRetain = NumRetain + 1
     else
        NumExportPE(DestinationPE(i)) = NumExportPE(DestinationPE(i)) + 1
     endif

  enddo

  ! Exchange data with other processes
  call mpi_barrier(mpi_comm_world,ierr)
  write(logunit,*)'   -[mpi] exchanging export lists with other processes'
  call flush(logunit)
  call mpi_allgather(NumExportPE, NumPEs, mpi_integer, ExchangeMatrix, NumPEs, mpi_integer, mpi_comm_world, ierr)

  ! Deconstruct ExchangeMatrix into NumImport - be wary that MyPE runs [0:n-1]
  do i=1,NumPEs
     NumImportPE(i) = ExchangeMatrix((i-1)*NumPEs + (MyPE+1))
  enddo
  NumImport = sum(NumImportPE)
  NumExport = sum(NumExportPE)
  write(logunit,*)'   -retaining:',NumRetain
  write(logunit,*)'   -exporting:',NumExport
  write(logunit,*)'   -importing:',NumImport
  call flush(logunit)

  ! Construct the import and export offset arrays
  do i=2,NumPEs
     ImportOffset(i) = ImportOffset(i-1) + NumImportPE(i-1)
     ExportOffset(i) = ExportOffset(i-1) + NumExportPE(i-1)
  enddo


end subroutine build_exchange_lists

! *************************************************

subroutine exchange_particles(itype)
  use logs
  use my_mpi
  use particles
  use exchange
  use smoothing_variables
  use header , only : BoxSize
  
  implicit none
  integer*4 :: itype, ivar, iPE
  integer*4 :: i, j, islab, ierr
  integer*4 :: nbytes, position, recv
  integer*4 :: level, nsend_levels, sendTask, recvTask
  integer*4 :: ii, jj, nsend, nrecv, sendTag,recvTag
  integer*4 :: NumPart_MyPE_Old, count
  integer*4, dimension(mpi_status_size) :: status

  write(logunit,*)'  -initiating particle exchange',NumExport
  call flush(logunit)

  ! Determine required buffer size
  nvars = (3+1)                        ! For coordinates, mass  
  nvars = nvars + 1                    ! For smoothing lengths
#ifdef FORCE_EOS
  if(itype == 0) nvars = nvars + 1     ! For equation of state flag
#endif
  nvars = nvars + nvar_by_type(itype)  ! For smoothing variables
  nbytes = nvars*4

  ! Construct the exchange buffer
  call allocate_export_buffer(itype,nbytes)
  
  write(mystring,'(D7.2)')dble(size(ExportBuf))/(1024.0d0*1024.0d0)
  write(logunit,*)'   -allocated '//trim(adjustl(mystring))//'Mb for export buffer'
  call flush(logunit)

  ! Reset this counter so we can compute the offset in ExportBuf
  NumExportPE(:) = 0

  do i=1,NumPart_ByPE(MyPE)
     islab = DestinationPE(i) 
     if(islab > 0 .and. islab <= NumPEs .and. islab /= (MyPE+1)) then
        position = nbytes*(ExportOffset(islab) + NumExportPE(islab)) 
        NumExportPE(islab) = NumExportPE(islab) + 1 
        call mpi_pack(readCoordinates(:,i),3,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
        call mpi_pack(readMass(i),1,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
        call mpi_pack(readSmoothingLength(i),1,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
#ifdef FORCE_EOS
        if(itype == 0) call mpi_pack(readEoS(i),1,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
#endif
        do ivar=1,nvar_by_type(itype)
             call mpi_pack(readSmoothVar(ivar,i),1,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
        enddo
     endif
  enddo
  
  ! Allocate import buffer
  call allocate_import_buffer(itype,nbytes)

  write(mystring,'(D7.2)')dble(size(ImportBuf))/(1024.0d0*1024.0d0)
  write(logunit,*)'   -allocated '//trim(adjustl(mystring))//'Mb for import buffer'
  call flush(logunit)
  
  ! Determine the number of possible communication pairs
  nsend_levels = 0
  do while(2**nsend_levels < NumPEs)
     nsend_levels = nsend_levels + 1
  enddo
  nsend_levels = 2**nsend_levels
  
  ! Implement bitwise XOR to determine send/receive pairs
  write(logunit,*)'   -[mpi] exchanging particle data with other processes'
  call flush(logunit)
  do level=0,nsend_levels-1
     sendTask = MyPE
     recvTask = IEOR(sendTask,level)
     sendTag = sendTask
     recvTag = recvTask
     if(recvTask < NumPEs)then
        ii = (nbytes*ExportOffset(recvTask+1)) + 1
        nsend = nbytes*NumExportPE(recvTask+1) 
        jj = (nbytes*ImportOffset(recvTask+1)) + 1
        nrecv = nbytes*NumImportPE(recvTask+1) 
        if(nsend == 0) ii = 1 ! Prevents seg fault when final processor isn't exporting
        if(nrecv == 0) jj = 1 ! Prevents seg fault when final processor isn't importing
        call mpi_sendrecv(ExportBuf(ii:ii),nsend,mpi_character,recvTask,sendTag,&
             ImportBuf(jj:jj),nrecv,mpi_character,recvTask,recvTag,mpi_comm_world,status,ierr)
     endif
  enddo

  write(logunit,*)'   -done import'
  call flush(logunit)

  ! Can now deallocate export buffer
  call deallocate_export_buffer()

  ! Allocate memory for this process' particles (both those read in and imported)
  NumPart_MyPE_Old   = NumPart_ByPE(MyPE)
  NumPart_ByPE(MyPE) = NumRetain + NumImport
  write(logunit,*)'   -new NumPart_MyPE =',NumPart_ByPE(MyPE)
  call allocate_particles(itype)

  ! Move non-exported local particles into new arrays
  count = 0
  do i=1,NumPart_MyPE_Old
     if(DestinationPE(i)==(MyPE+1))then
        count = count + 1
        Coordinates(:,count)   = readCoordinates(:,i)
        Mass(count)            = readMass(i)
        SmoothingLength(count) = readSmoothingLength(i)
#ifdef FORCE_EOS
        if(itype==0) EoS(count) = readEoS(i)
#endif
        do ivar=1,nvar_by_type(itype) 
            SmoothVar(ivar,count) = readSmoothVar(ivar,i)
        enddo
     endif
  enddo

  ! Now free up the read_particles array
  call deallocate_read_particles()

  ! Now unpack the import buffer into the remainder of the array
  position = 0
  call flush(logunit)
  do i=1,NumPEs
     if(NumImportPE(i) > 0)then
        do j=1,NumImportPE(i)
           count = count + 1
           call mpi_unpack(ImportBuf,size(ImportBuf),position,Coordinates(:,count),3,mpi_real,mpi_comm_world,ierr)        
           call mpi_unpack(ImportBuf,size(ImportBuf),position,Mass(count),1,mpi_real,mpi_comm_world,ierr)        
           call mpi_unpack(ImportBuf,size(ImportBuf),position,SmoothingLength(count),1,mpi_real,mpi_comm_world,ierr)        
#ifdef FORCE_EOS
           if(itype==0) call mpi_unpack(ImportBuf,size(ImportBuf),position,EoS(count),1,mpi_real,mpi_comm_world,ierr)        
#endif
           do ivar=1,nvar_by_type(itype) 
              call mpi_unpack(ImportBuf,size(ImportBuf),position,SmoothVar(ivar,count),1,mpi_real,mpi_comm_world,ierr)        
           enddo
        enddo
     endif
  enddo
  
  ! Can now also free up some memory
  call mpi_barrier(mpi_comm_world,ierr)
  call deallocate_import_buffer()
  call deallocate_destination_array()
  call deallocate_exchange_counters()

end subroutine exchange_particles


! *************************************************

subroutine exchange_ghosts(itype)
  use logs
  use my_mpi
  use parameter_module
  use camera_module, only : z0, L
  use exchange
  use grid
  use particles
  use smoothing_variables
  use header, only : BoxSize
  implicit none
  integer*4 :: itype, ThisPE, ipart, i, j, ierr, ivar, iPE
  integer*4 :: Noff, nbytes
  integer*4 :: level, nsend_levels, sendTask, recvTask
  integer*4 :: sendTag, recvTag
  integer*4 :: count, nsend, nrecv, recv, ii, jj, position
  integer*4, dimension(mpi_status_size) :: status
  real*4 :: z, h, zmin, zmax
  real*4    :: hmin

  write(logunit,*)'   -building list of ghosts to be exported'
  call flush(logunit)

  ! We'll need these arrays again
  call allocate_exchange_counters()
  call allocate_destination_array()
  NumRetain = 0
  NumImportPE(:) = 0
  NumExportPE(:) = 0
  ExchangeMatrix(:) = 0
  ImportOffset(1) = 0
  ExportOffset(1) = 0

  ! Define a minimum smoothing length
#ifdef CUBE
  hmin = sqrt(3.0/8.0)*pixelsize
#else
  hmin = sqrt(1.0/2.0)*pixelsize
#endif

  ! Determine required size of export buffer
  do ThisPE = 0,NumPEs-1

     ! Don't want to export to self!
     if(ThisPE /= MyPE)then

        Noff = (NZMaster + (ThisPE-1)*NZLocal)
#ifdef CUBE
        GhostExtent(1) = GridExtent(1,3)  + real(Noff)*CellSize(3) 
        GhostExtent(2) = GhostExtent(1)   + real(NZLocal)*CellSize(3)
#else
        GhostExtent(1) = GridExtent(1,3)  + real(Noff)*CellSize(2) 
        GhostExtent(2) = GhostExtent(1)   + real(NZLocal)*CellSize(2)
#endif
        do ipart=1,NumPart_ByPE(MyPE)
           z = Coordinates(3,ipart) 
           h = max(SmoothingLength(ipart),hmin)
           zmax = z+h
           zmin = z-h
           if(zmax > GhostExtent(1) .and. zmin < GhostExtent(2)) then
              ! Must export this particle to this process
              NumExportPE(ThisPE+1) = NumExportPE(ThisPE+1) + 1
           endif
        enddo

     endif
  enddo

  ! Exchange export data with other processes
  write(logunit,*)'   -[mpi] exchanging ghost export lists with other processes'
  call flush(logunit)
  call mpi_allgather(NumExportPE, NumPEs, mpi_integer, ExchangeMatrix, NumPEs, mpi_integer, mpi_comm_world, ierr)

  ! Deconstruct ExchangeMatrix into NumImport - be wary that MyPE runs [0:n-1]
  do i=1,NumPEs
     NumImportPE(i) = ExchangeMatrix((i-1)*NumPEs + (MyPE+1))
  enddo

  NumImport = sum(NumImportPE)
  NumExport = sum(NumExportPE)
  write(logunit,*)'   -exporting:',NumExport
  write(logunit,*)'   -importing:',NumImport
  call flush(logunit)
  call flush(logunit)

  ! Construct the import and export offset arrays
  do i=2,NumPEs
     ImportOffset(i) = ImportOffset(i-1) + NumImportPE(i-1)
     ExportOffset(i) = ExportOffset(i-1) + NumExportPE(i-1)
  enddo

  ! Set how many bytes to offset by
  nvars = (3+1)                     ! For coordinates, mass
  nvars = nvars + 1                 ! For smoothing lengths
#ifdef FORCE_EOS
  if(itype==0) nvars = nvars + 1    ! For equation of state flag
#endif
  nvars = nvars + nvar_by_type(itype)
  nbytes = nvars*4

  ! Allocate the export buffer
  call allocate_export_buffer(itype,nbytes)  

  write(mystring,'(F7.2)')real(size(ExportBuf))/(1024.0*1024.0)
  write(logunit,*)'   -allocated '//trim(adjustl(mystring))//'Mb for ghost export buffer'
  call flush(logunit)

  ! Now populate the export buffer
  position = 0
  do ThisPE = 0,NumPEs-1

     ! Don't want to export to self!
     if(ThisPE /= MyPE)then

        Noff = (NZMaster + (ThisPE-1)*NZLocal)
#ifdef CUBE
        GhostExtent(1) = GridExtent(1,3)  + real(Noff)*CellSize(3) 
        GhostExtent(2) = GhostExtent(1)   + real(NZLocal)*CellSize(3)
#else
        GhostExtent(1) = GridExtent(1,3)  + real(Noff)*CellSize(2) 
        GhostExtent(2) = GhostExtent(1)   + real(NZLocal)*CellSize(2)
#endif
        do ipart=1,NumPart_ByPE(MyPE)
           z = Coordinates(3,ipart) 
           h = max(SmoothingLength(ipart),hmin)
           zmax = z+h
           zmin = z-h
           if(zmax > GhostExtent(1) .and. zmin < GhostExtent(2)) then

              call mpi_pack(Coordinates(:,ipart),3,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
              call mpi_pack(Mass(ipart),1,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
              call mpi_pack(SmoothingLength(ipart),1,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
#ifdef FORCE_EOS
              if(itype==0) call mpi_pack(EoS(ipart),1,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
#endif
              do ivar=1,nvar_by_type(itype) 
                   call mpi_pack(SmoothVar(ivar,ipart),1,mpi_real,ExportBuf,size(ExportBuf),position,mpi_comm_world,ierr)
              enddo

           endif
        enddo
     endif
  enddo


  ! Allocate import buffer
  call allocate_import_buffer(itype,nbytes)

  write(mystring,'(F7.2)')real(size(ImportBuf))/(1024.0*1024.0)
  write(logunit,*)'   -allocated '//trim(adjustl(mystring))//'Mb for ghost import buffer'
  call flush(logunit)

  ! Determine the number of possible communication pairs
  nsend_levels = 0
  do while(2**nsend_levels < NumPEs)
     nsend_levels = nsend_levels + 1
  enddo
  nsend_levels = 2**nsend_levels
       
  ! Implement bitwise XOR to determine send/receive pairs
  write(logunit,*)'   -[mpi] exchanging particle data with other processes'
  call flush(logunit)
  do level=0,nsend_levels-1
     sendTask = MyPE
     recvTask = IEOR(sendTask,level)
     sendTag = sendTask
     recvTag = recvTask
     if(recvTask < NumPEs)then
        ii = (nbytes*ExportOffset(recvTask+1)) + 1
        nsend = nbytes*NumExportPE(recvTask+1) 
        jj = (nbytes*ImportOffset(recvTask+1)) + 1
        nrecv = nbytes*NumImportPE(recvTask+1) 
        if(nsend == 0) ii = 1 ! Prevents seg fault when final processor isn't exporting
        if(nrecv == 0) jj = 1 ! Prevents seg fault when final processor isn't importing
        call mpi_sendrecv(ExportBuf(ii:ii),nsend,mpi_character,recvTask,sendTag,&
             ImportBuf(jj:jj),nrecv,mpi_character,recvTask,recvTag,mpi_comm_world,status,ierr)
     endif
  enddo

  ! Can now deallocate export buffer
  call deallocate_export_buffer()

  ! Allocate memory for this process' imported ghost particles
  NumGhost_ByPE(MyPE) = NumImport

  ! Broadcast local particle number to each processor
  call mpi_barrier(mpi_comm_world,ierr) 
  do iPE=0,NumPEs-1
     call mpi_bcast(NumGhost_ByPE(iPE),1,mpi_integer,iPE,mpi_comm_world,ierr)
  enddo

  call allocate_ghost_particles(itype)

  ! Now unpack the import buffer into the remainder of the array
  write(logunit,*)'   -unpacking received particle data'
  call flush(logunit)
  position = 0
  count = 0
  do i=1,NumPEs
     if(NumImportPE(i) > 0)then
        do j=1,NumImportPE(i)
           count = count + 1
           call mpi_unpack(ImportBuf,size(ImportBuf),position,Ghost_Coordinates(1:3,count),3,mpi_real,mpi_comm_world,ierr)        
           call mpi_unpack(ImportBuf,size(ImportBuf),position,Ghost_Mass(count),1,mpi_real,mpi_comm_world,ierr)        
           call mpi_unpack(ImportBuf,size(ImportBuf),position,Ghost_SmoothingLength(count),1,mpi_real,mpi_comm_world,ierr)        
#ifdef FORCE_EOS
           if(itype==0) call mpi_unpack(ImportBuf,size(ImportBuf),position,Ghost_EoS(count),1,mpi_real,mpi_comm_world,ierr)  
#endif
           do ivar=1,nvar_by_type(itype) 
                 call mpi_unpack(ImportBuf,size(ImportBuf),position,Ghost_SmoothVar(ivar,count),1,mpi_real,mpi_comm_world,ierr)
           enddo
        enddo
     endif
  enddo
  write(logunit,*)'    -done'
  call flush(logunit)

  ! Can now also free up some memory
  call mpi_barrier(mpi_comm_world,ierr)
  call deallocate_import_buffer()
  call deallocate_exchange_counters()

end subroutine exchange_ghosts
